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We describe a method for imaging 3D objects in a tomographic configuration implemented 
by training an artificial neural network to reproduce the complex amplitude of the exper¬ 
imentally measured scattered light. The network is designed such that the voxel values of 
the refractive index of the 3D object are the variables that are adapted during the training 
process. We demonstrate the method experimentally by forming images of the 3D refractive 
index distribution of cells. 

The learning approach to imaging we describe in this paper is related to adaptive techniques 
in phased antenna arrays 1 and inverse scattering 2 3 . In the optical domain an iterative approach was 
demonstrated by the Sentenac groupie who used the coupled dipole approximation^ for modelling 
light propagation in inhomogeneous media (a very accurate method but computationally intensive) 
to simulate light scattering from small objects (1 pm x 0.5 pm) in a point scanning microscope 
configuration. Very recently an iterative optimization method was demonstrated 7 * for imaging 3D 
objects using incoherent illumination. Our method relies on digital holography 8 9 to record the 
complex amplitude of the field. We use the Beam Propagation Method (BPM) 10 11 to model the 
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scattering process and the error back propagation method 12 to train the system. At the end of the 
training process the network discovers a 3D index distribution that is consistent with the experi¬ 
mental observations. We experimentally demonstrate the technique by imaging polystyrene beads 
and HeLa and hTERT-RPEl cells. 

Experimental Setup 

A schematic diagram of the experimental setup is shown in Figure 1. It is a holographic tomog¬ 
raphy system 13 , in which the sample is illuminated with multiple angles and the scattered light 
is holographically recorded. Several variation of the holographic tomography system have been 
demonstrated before 14 17 . The optical arrangement we used is most similar to the one described 
by Choi et al. u . The samples to be measured were prepared by placing polystyrene beads and 
cells between two glass cover slides. The samples were illuminated with a continuous collimated 
wave at 561 ram at 80 different angles. The amplitude and phase of the light transmitted through 



Figure 1: Experimental setup (BS: Beam Splitter, GM: Galva Mirror, L: Lens, OB: Objective, M: 
Mirror) 
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the sample was imaged onto a 2D detector where it was holographically recorded by introducing 
a reference beam. These recordings constitute the training set with which we train the computa¬ 
tional model that simulates the experimental setup. We construct the network using the BPM. The 
inhomogeneous medium (beads or cells) is divided into thin slices along the propagation direction 
(z). The propagation through each slice is calculated as a phase modulation due to the local trans¬ 
verse index variation followed by propagation in a thin slice of a homogenous medium having the 
average value of the index of refraction of the sample. 

Methodology 

A schematic description of the BPM simulation is shown in Figure 2. The straight lines connecting 
any two circles represent multiplication of the output of the unit located in the /-th layer of the net¬ 
work at x = ni8,y = miS by the discretized Fresnel diffraction kernel e^ n ^~ n ^+ 1 ' )52+ ^~ m ^+ 1 ' >s2 ^ XSz 
where «/ and m/ are integers and A is the wavelength of light. 8 is the sampling interval in the trans¬ 
verse coordinates (x, y ) whereas 8 Z is the sampling interval along the propagation direction 2 . The 
circles in the diagram of Figure 2 perform a summation of the complex amplitude of the signals 
converging to each circle and also multiplication of this sum by e i( 27TAnS ^ z )/ x _ A n(x, y, z ) is the 
unknown 3D index perturbation of the object. 

In the experiments the network has 420 layers with A n(x, y, z) being the adaptable variable. 
In contrast with a conventional neural network, the output of the layered structure in Figure 2 is 
a linear function of the input complex field amplitude. However, the dependence of the output is 
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Figure 2: Schematic diagram of object reconstruction by learning the 3D index distribution that 
minimizes the error e, defined at the mean squared difference between the experimental measure¬ 
ment and the prediction of a computational model based on the beam propagation method (BPM). 


nonlinearly related to A n(x, y, z). The BPM can be trained using steepest descent exactly as the 
back propagation algorithm in neural networks 18 20 . Specifically, the learning algorithm carries out 
the following minimization: 

mini- 1 - \\E k (An) - M k (An)\\ 2 + TS(An)} 

An k =1 


subject to 0 < An 

In the above expression E k {An) is the current prediction of the BPM network for the output when 
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the system is illuminated with the fc-th beam and M k (An) is the actual measurement obtained by 
the optical system. A n indicates the estimate for the index perturbation due to the object. The 
term S(Ah) is a sparsity constraint 21 to enhance the contrast while r is a parameter that can be 
tuned to maximize image quality by systematic search. The positivity constraint takes advantage 
of the assumption that the index perturbation is real and positive. The optimization is carried out 
iteratively by taking the derivative of the error with respect to each of the adaptable parameters 
following steepest descent 


An —>• An — 


a y, de^ OSjAn) 
fc ^ <9An T <9 A n 


where e k = \\E k (An) — M k (An) || is the error, a is a constant and the change in An is proportional 
to the error and its derivative. This is achieved efficiently via a recursive computation of the 
gradient, which is the back propagation part of our learning algorithm. 


Results 

We first tested the system with polystyrene beads encapsulated between two glass slides in im¬ 
mersion oil. The sample was inserted in the optical system of Figure 1 and 80 holograms were 
recorded by illuminating the sample at 80 distinct angles uniformly distributed in the range -45 
degrees to +45 degrees. The collected data is the training set for the 420-layer BPM network 
which simulates a physical propagation distance of 30 (xm and transverse window 37/jm x 37 fim 
(5 X = <) y = 72 nm). The network was initialized with the standard filtered back projection recon¬ 
struction algorithm (Radon transform) 22 and the resulting 3D images before and after 100 iterations 
are shown in Figure 3. The final image produced by the learning algorithm is an accurate repro- 
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Figure 3: Reconstruction of two 10/im beads of refractive index 1.588 at A = 561nm in im¬ 
mersion oil with n 0 = 1.516. (a)-(c) x-y, y-z and x-z slices using the inverse Radon transform 
reconstruction, (d)-(f) the same slices for our learning based reconstruction method. 

duction of the bead shape. 

A sample of a HeLa cell was also prepared and the same procedure was followed to obtain 
a 3D image. The results are shown in Figure 4 where the error function is plotted as a function 
of iteration number. In this instance, the system was initialized with a constant but nonzero value 
(An = 0.007). Also shown in Figure 4 are the results obtained when the system was initialized 
with the Radon reconstruction from the same data. After 100 iterations both runs yield essentially 
identical results. Notice that the error in the final image (after 100 iterations) is significantly lower 
than to the error of the Radon reconstruction. This is also evident by visual inspection of the 
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images in Figure 4 where the artifacts due to the missing cone 23 and diffraction 14 are removed by 
the learning process. 

In general, optical 3D imaging techniques rely on the assumption that the object being im¬ 
aged does not significantly distort the illuminating beam. This is assumed for example in Radon or 
diffraction/holographic tomography. In other words, these 3D reconstruction methods rely on the 
assumption that the measured scattered light consists of photons that have only been scattered once 
before they reach the detector. The BPM, on the other hand, allows for multiple forward scattering 
events. The only simplification is that reflections are not taken into account; these could eventually 
be incorporated in the network equation without fundamentally altering the approach described in 
this paper. Since biological tissue is generally forward scattering, BPM can be a good candidate 
to model propagation of thick biological samples and this may be the most significant advantage 
of the learning approach. To demonstrate this point, we prepared two glass slides with a random 
distribution of hTERT-RPEl cells (immortalized epithelial cells from retina) on each slide. When 
we attach the two slides together, we can find locations where two cells are aligned in 2 , one on 
top of the other. Figure 5 (a)-(e) shows the image of such a stack of two cells produced with a 
direct inversion using the Radon transform. Figure 5 (f)-(j) shows the same object imaged with 
the proposed learning algorithm. The learning method was able to distinguish the two cells where 
the Radon reconstruction merged the two into a single pattern due to the blurring in z which is a 
consequence of the missing cone. 

In conclusion, we have demonstrated a neural-network based algorithm to solve the optical 
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Figure 4: Comparison of the proposed method initialized by inverse Radon transform (left) versus 
initialization with a constant value (An = 0.007) (right), (a),(e) pixel error energy fall-off for 
initialized with inverse Radon and constant value, respectively. The horizontal doted line shows 
the inverse Radon performance for comparison, (b)-(d), x-y, y-z and x-z stacks for respectively the 
first, tenth and hundredth iteration of the proposed method initialized by inverse Radon, (d)-(f), 
the same figures for the proposed method initialized by constant value. 
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Figure 5: Images of two hTERT-RPEl cells, x-y slices corresponding to different depths of re¬ 
spectively +9, +6, +3, 0 and -3 microns (positive being toward the detector) from the focal plane of 
the lens OB2 in Figure 1 for: (a)-(e) the inverse Radon transform based reconstruction and (f)-(j) 
the same slices for our learning based reconstruction method. 
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phase tomography problem and have applied it to biological (HeLa and hTERT-RPEl cells) and 
synthetic (polystyrene beads) samples. The experimental measurements were performed with a 
conventional collimated illumination phase tomography setup, with coherent light and holograms 
were recorded off-axis. The sample scattering potential was modeled as a neural network imple¬ 
menting a forward beam propagation method. The network is organised in neurones layers, each 
one of them representing an x-y plane in the BPM. The output of the network is compared to the 
experimental measurements and the error is used to correct the weights (representing the refractive 
index contrast) in the neurones using standard error back propagation techniques. The algorithm 
yields images of better quality than tomographic reconstruction (Radon). In particular, the missing 
cone artefact is efficiently removed, as well as parasitic granular structure. We have shown that 
whether starting from a constant initial guess for the refractive index or with a conventional Radon 
tomographic image, the method essentially converges to the same result after 100 iterations. This 
approach opens rich perspectives for active correction of scattering in biological sample; in partic¬ 
ular, it has the potential of increasing the resolution and the contrast in fluorescent and two-photon 
imaging. 
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